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Abstract 

The paper describes a physical model and numerical algorithm for modeling Type la supernova 
(SNIa) explosions in three-dimensions and presents first results of modeling a defiagration explosion 
in a non-rotating, Chandrasekhar-mass carbon-oxygen (CO) white dwarf. Simulations show that 
the turbulent fiame speed grows exponentially, reaches ~ 30% of the speed of sound, and then 
declines as the large-scale turbulence is frozen by expansion. The freezing of turbulent motions 
appears to be a crucial physical mechanism regulating the rate of deflagration in SNIa. The energy 
of the explosion is comparable to that of a typical SNIa. However, the presence of the outer layer 
of unburned CO and the formation of intermediate mass elements and pockets of unburned CO 
near the center pose problems for the modeling of SNIa spectra. Delayed detonation is a way to 
alleviate these problems and to produce consistent spectra. 



1. Introduction 



Type la Supernovae (SNIa) are thought to be thermonuclear explosions of Chandrasekhar-mass 
carbon-oxygen white dwarfs (CO-WD). They play an important role in nucleosynthesis (Nomoto 
et al. 1984) and the balance of energy in interstellar medium, and they are important distance 
indicators used in cosmology (Perlmuter et al., 1999; Riess et al. 2000). Our understanding 
of SNIa explosions is far from complete. To predict the explosion outcome, one needs to know 
how thermonuclear burning propagates inside the exploding star. The mechanisms and the speed 
of thermonuclear burning in SNIa remain an unsolved theoretical problem. A recent review by 
Hillebrandt & Niemeyer (2000) contains a large list of relevant publications. 

Burning can propagate in a supernova either as a flame caused by the electron heat conduction 
or as a detonation in which reactions are triggered by a shock. Both the detonation speed and 
the laminar flame speed Si in SNIa are known. However, in the gravitational field of the star, 
the laminar flame becomes unstable with respect to the Rayleigh- Taylor (RT) instability (Nomoto 
et al. 1976). This leads to a turbulent deflagration with some unknown turbulent flame speed 
St > Si. 

Observations indicate that burning of a Chandrasekhar-mass CO-WD must start as a defla- 
gration. Spectra of SNIa at maximum light show the presence of large amounts of intermediate 
mass elements such as Si, S, Mg (Pskovskii 1977, Branch 1981). These elements can only be 
synthesized at densities less than ~ lO^g/cc, where the temperature of burned products stays 
below ~ (4 — 5) X lO^K. When burning starts as a deflagration, the initially high-density WD 
expands while it burns, its density decreases, and this creates the necessary conditions for the 
production of intermediate mass elements. In one-dimensional deflagration models, St is a free 
parameter (Nomoto et al. 1976, Woosley & Weaver 1986). Delayed detonation models assiime 
that a deflagration undergoes a transitions to a detonation at low densities (Khokhlov 1991). In 
these models, both St and the moment of a deflagration-to-detonation transition (DDT) are free 
parameters. To fit observations, burning must be initially slow in order to expand the WD signif- 
icantly. Then burning must accelerate quickly to incinerate expanding outer layers rapidly. This 
requires very carefully tuning St as a function of time in deflagration models. Delayed-detonation 
models accomplish the incineration by switching to a supersonic detonation mode of burning. 

Two-dimensional simulations of SNIa explosions show that burning leads to a significant ex- 
pansion of a WD. However, the deflagration speed found in all of the simulations is too small to 
lead to a powerful explosion (Livne, 1993; Arnett Sz Livne, 1994; Reinecke et al., 2000). Whether 
this is an inherent difficulty of the defiagration model or simply a deficiency of two-dimensional 
simulations is a question. Two-dimensional simulations lack a key physical ingredient of turbu- 
lent burning - turbulent energy cascade that is directed in three dimensions from large to small 
spatial scales. Simulations of burning in a vertical column in a uniform gravitational field show 
that the turbulent flame speed grows faster in three dimensions and eventually reaches a turbulent 
steady state where St becomes independent of Si and of the details of burning on small scales 
(Khokhlov, 1995; hereafter K95). In two dimensions, St grows more slowly and depends on the 
value of St. However, turbulence in the exploding star is affected by expansion on spatial scales 
C where the characteristic RT timescale tj^t — (5^)^^^ is comparable or exceeds the expansion 
timescale Tg R/Ug, where R is stellar radius and Ug is the (time-dependent) expansion velocity. 
Expansion tends to freeze large-scale turbulent motions and thus decreases the turbulent flame 
speed (K95). The actual St in a supernova must depend on the competition of the RT instability, 
turbulent energy cascade, and expansion. To predict St, these effects must be modeled in three 
dimensions. This paper describes a physical model and numerical algorithm for modeling SNIa 
explosions in three-dimensions and presents flrst results of modeling a deflagration explosion in a 
nonrotating, Chandrasekhar-mass CO-WD. 
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2. Physical and numerical model 

2.1. Fluid dynamics 
Stellar matter is described by the Euler equation for an inviscid fluid 

V • (pU) , 

V-(pUU)-VP + pg, (1) 
V-(U {E + P))+p\J-g + pq, 

where p, E = Ei + pU'^/2, Ei, U, g and q are the mass density, energy density, internal energy 
density, velocity of matter, gravitational acceleration, and nuclear energy release rate per unit mass, 
respectively. The acceleration of gravity is computed using an angle-averaged density distribution. 
The equation of state (EOS) of degenerate matter includes contributions from ideal Fermi-Dirac 
electrons and positrons, equilibrium Planck radiation, and ideal ions. Pressure P = P{p, Ei,Ye, Yi) 
and temperature T = T{p, Ei,Yf.,Yi) are determined by the EOS as functions of p, E^, electron 
mole fraction Yg, and the mean mole fraction of ions Yi. 



dp 

di 
dpU 

dE 
'dt 



2.2. Flame propagation 



The flame is advanced using a flame-capturing technique which mimics a flame with a pre- 
scribed normal speed (K95). A scalar variable /, such that / = in the unburned matter and 
/ = 1 in the material which has passed through the flame, obeys a reaction-diffusion equation 

^ + U-V/ = irvV + i?, (2) 

with artificial reaction and diffusion coefficients 



R 



C = const., if /o < / < 1; 

0, otherwise, (3) 

K = const. 



with some < /o < 1. In this paper, /o = 0.3. Equation (2) has a solution f{x ~ St) which 
describes a reaction front that propagates with the speed S = {KC / fo)^^'^ and has a thickness 
5 ~ [K/Cf/^ (see Appendix of (K95)). If K and C in (3) are set to 

K = S{f5Ar)^o, C'=^/^, (4) 

with (3 = constant, the front propagates with a prescribed speed S and spreads onto several 
computational cells of size Ar. The choice of /? = 1.5 spreads the flame on ~ 3 — 4 cells. Making the 
front narrower is not practical, since the fluid dynamics algorithm spreads contact discontinuities 
on ~ 4 cells. The energy-release rate inside the front is defined as 

^=<ifJt^ (5) 
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where qf, defined below in Section 3.4, is the nuclear energy release inside the front. Due to the 
energy release, both density and velocity vary across the front, and fluid motions are generated 
which advect the flame relative to the computational mesh. The flame-capturing technique de- 
scribed above advances the reaction front relative to the fuel with the speed practically independent 
of the orientation of the front on the mesh, fluid motions, and resolution. 

2.3. Burning on small scales 

A turbulent flame in a vertical column of width C subjected to an acceleration of gravity g 
propagates in a steady-state with the speed St — {).b\J AgC independent of the laminar speed S*/, 
where A = (po — Pi)/{po + pi) is the Atwood number, and po, pi are the densities ahead and 
behind the flame front, respectively (K95). The time to reach a steady state is ~ 2C/St, and the 
thickness of the flame brush is ^ 2C. To describe burning on scales that are not resolved in the 
simulations, the value of S in (4) is taken as 

S = max{Si,Ssub) , (6) 

with 

Ssub = 0.5a/A aAr (si |g| + S2 max (0, n • g)) (7) 

based on the assumption that burning on small scales is driven by the RT instability at these 
scales and that at C « Rwd burning adjusts to a local steady state (Livne 1993, Arnctt & Livne 
1994, K95). In (7), a ~ 1 determines the driving scale, C ~ aAr, in terms of the computational 
cell size Ar, and n = Vf/f is the unit vector normal to the flame front and directed towards 
the products of burning. The choice of si = 1, S2 = gives a subgrid flame speed that is 
independent of the orientation of the flame front relative to the direction of gravity. The choice 
si = 0, S2 = 1 gives a directional-dependent speed that reduces to Si when the flame moves in 
the direction of gravity. Finally, si = S2 = gives the laminar flame speed. In this paper, only 
the limiting, angle-independent and laminar prescriptions for S were used. The Atwood number 
A ~ (7 — 1) pgj / 27P was calculated from conservation of enthalpy, assuming constant pressure 
across the flame, where 7 and p = {pi + po)/2 are average values of the adiabatic index and mass 
density. Under the conditions of interest here, A — 0.4. The laminar flame speed Si, taken 

according to Timmes & Woosley (1992) and Khokhlov et al. (1997), is Si ~ (10"^ - 10"^) x 
for p ~ 10^ — 10''g/cc, where ~ 5 x 10*cm/s is the sound speed. 

2.4- Nuclear kinetics 

There are three distinct stages of carbon burning. First, the ^^C -|- ^^C reaction leads to 
the consumption of ^^C and formation of mostly "^^Ne, '^^Mg, protons and a-particles. Then the 
onset of the nuclear statistical quasi-equilibrium (NSQE) takes place, during which Si-group (in- 
termediate mass) elements are formed. Finally, Si-group elements are converted into the Fe-group 
elements and the nuclear statistical equilibrium (NSE) sets in. The NSQE and NSE relaxation 
involves hundreds of species from carbon to zinc and thousands of reactions of these nuclei with 
protons, neutrons, and a-particles. The timescales of all these stages strongly depend on temper- 
ature. The timcscalc of carbon consumption is much shorter than the quasi-equilibrium timescale 
Tnsqe which, in turn, is much shorter than the NSE timescale Tnse- 

A four-equation kinetic scheme is used here to describe the energy release, synthesis of Si- 
and Fe-group elements, and the neutronization of NSE matter. A similar scheme has been used in 
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one-dimensional simulations of SNIa (Khokhlov 1991). The kinetic equation for the mole fraction 
of carbon Yc 

^ = -pA{T,) eM-Q/TlhY^ , (8) 

describes carbon consumption through the major reactions ^^C {^"^C, p)'^^Na{p,j)'^'^Mg and 
i2^(^i2(^^ ^Hey^Ne with the branching ratio 1. Here Q = 84.165, Tg^ = Tg/fl + O.OeTTg), 
where A(Tg) is a known function (Fowler et al. 1975). 

Most of the nuclear energy is released during the carbon exhaustion stage and the subsequent 
synthesis of the Si-group nuclei. The energy release or consumption due to the transition from Si- 
group to Fe-group nuclei (NSE relaxation) is less than 10%. Therefore, the nuclear energy release 
rate is approximated as 

dqn „ dVc , qnse — Qn /r>\ 

^=-^^^+^-^' 

where Qnit) is the binding energy of nuclei per unit mass, Qc = 4.48 X IQ-*"^ ergs g~^ mol~^ is the 
energy release due to carbon burning (8), and qnse{p,T,Yf,) is the binding energy of matter in the 
state of NSE. 

The equation 

dSnse 1 ^nse (10) 

dt ^nse 

traces the onset of NSE, where 5nse = in the unburned matter and Snse = 1 in the NSE products. 
Intermediate mass elements are expected where Yc{t = oo) ~ and 5nse(i = oo) < 1. Fe-group 
elements are expected where Yc{t = oo) c± and Snseit = oo) c± 1. The 'e-folding' NSQE and NSE 
timescales 

T„,,e = exp (149.7/r9 - 39.15) s , 

Tnse = exp (179.7/rg - 40.5) s , ^ ^ 

approximate the results of the detailed calculations of carbon burning. The original kinetic scheme 
(Khokhlov 1991) also contained the equation for the mean ion mole fraction YJ. To simplify the 
scheme in this paper, variations of Yi due to nuclear reactions were neglected since ions make a 
small contribution to the equation of state. A constant Yi = 0.07 was used instead, which is an 
average of Y^ in the unburned and typical NSE matter in the SNIa explosion conditions. 
Neutronization is described by the equation for the electron mole fraction Y^, 

dY 

= -R^{p,T,Y,) . (12) 

The corresponding term describing neutrino energy losses, — g^(p, T, Ye), is added to the energy 
conservation equation in (1). Values of qnse, qw, and were computed assuming NSE distribution 
of individual nuclei, as described in (Khokhlov 1991). Recently, there has been an important 
development in theoretical computations which shows significantly smaller electron capture and 
/3-decay rates in stellar matter (Langanke et al. 1999 and references therein). Following Martinez- 
Pinedo et al. (2000), this effect was approximately taken into account by decreasing q^ and 
by a factor of 5. This should be sufficient to account for changes in q^se and in the corresponding 
nuclear energy release qnse — qn{t = 0) caused by variations of Ye- In future simulations, we plan 
to use q-uj and R^ based on the new individual capture rates. 

Despite its simplicity, the kinetic scheme reasonably well describes all major stages of carbon 
burning. In particular, it takes into account the important effect of energy release or absorption 
caused by changes in NSE composition when matter expands or contracts (these changes happen 
on a quasi-equilibrium rather than on the equilibrium timescale). When the WD expands, this 
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effect adds ~ 50% to the energy initially released by burning at high densities. Using even the 
simplest 13-species a-network to account for this effect would be prohibitively expensive in three- 
dimensional hydrodynamical simulations. 

To couple the kinetic scheme and the front tracking algorithm (Section 2.2), qf in (5) is taken 
Qf = Qnse — ?n(0) at densities p > 2 x 10'' g/cc, qf = QcYc{t = 0) at p < 5 x 10^ g/cc, and is 
linearly interpolated between the two values in the density range 5 x 10^ — 2 x 10^ g/cc. The carbon 
mole fraction Yc inside the front is decreased in proportion to the increase of /. Thus, at high 
densities, both carbon consumption and the NSQE relaxation take place inside the flame front. At 
low densities where the NSQE stage of burning is slow, NSQE relaxation takes place outside the 
flame front. 

2.5. Numerical method 

Fluid dynamics equations (1) for p, U, E, and advection parts of equations (2), (8)-(10) 
and 912) for chemical variables /, Yc, qn, Snse and Ye are solved using an explicit, second-order, 
Godunov-type, adaptive mesh refinement algorithm ALLA (Khokhlov, 1998). A Riemann solver 
for an arbitrary EOS (Colella & Glaz, 1985) is used to evaluate fluxes at cell interfaces. Left and 
right states are computed using a monotone piecewise-linear reconstruction (van Leer 1979). A 
small amount of an artificial diffusion flux (Lapidus, 1967) is added to Euler fluxes following Colella 
and Woodward (1984). This feature is used only in the vicinity of strong shocks. Multidimensions 
are treated using directional splitting. The diffusion term in (2) is treated using an explicit, 
second-order finite differencing. 

ALLA uses a fully threaded tree (FTT) adaptive mesh refinement to refine the computational 
mesh on the level of individual cells (Khokhlov 1998; Khokhlov & Chtchelkanova, 1999). FTT is a 
tree with inverted thread pointers directed from children to neighbors of parent cells. As opposed 
to conventional tree structures, all operations on the FTT, including tree modifications (refinement 
and unrefinement) are 0{N), where N is a number of cells. All of these operations can be performed 
in parallel. FTT has a 2| integers-per-cell memory overhead, compared with 8 integers for an 
ordinary tree or 14 integers for a tree that stores connectivity information. Computations using 
the FTT are virtually as fast (per computational cell) as those on a regular grid. A massively 
parallel version of the code uses a space-filling curve approach to maintain data locality, achieve 
load balance, and minimize interprocessor communications. ALLA has been extensively tested and 
used in various combustion problems involving shocks, flames, turbulence, and their interactions 
(Gamezo et al., 1999,2000; Khokhlov Sz Oran, 1999; and references therein) and in astrophysics 
(Khokhlov et al., 1999). 

Mesh refinement was used to refine shocks, contact discontinuities, and regions with large 
gradients of p, P, and /. The expressions for the corresponding shock, ^g, discontinuity, ^c, and 
gradient refinement indicators, ^p, ^p, (^f, are given by the equations (11), (12) and (13) of Khokhlov 
(1998), respectively. In addition, refinement was introduced for the gradients of tangential velocity 
using an indicator 

with a = 0.05. Refinement was initiated in places where the maximum of all indicators, 

^ = max(^s, Cp, Cp, C/, C«) , (14) 

was larger than the threshold value ^ > 0.5. The mesh was unrefined in places where ^ was less 
than the threshold value ^ < 0.05. The cell size Ar is related to the level I of the cell in the tree 
as Ar = L/2K The mesh can be characterized by the minimum and maximum levels of leaves 
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(unsplit cells) in the tree, Imin aiid Imax: which were predefined at the beginning of a computation. 
During a computation, mesh refinement was allowed for cells with Imin < I < Imax- 

Computations were carried out on a 128 250MHz RIOOOO SGI Origin 2000 at the Naval 
Research Laboratory (NRL). On this machine, ALLA advanced c± 2 x 10^ cells at the rate of one 
timestep per second per processor on pure fluid-dynamic, three-dimensional problems with 7 held 
constant. With the addition of a degenerate matter EOS, flame capturing, and nuclear kinetics, 
the speed decreased to ~ 6 x 10^ cell-steps/s/processor. 

2.6. Problem setup and initial conditions 

Most of the computation were performed for one octant of the WD assuming mirror symmetry 
along the x = 0, y = and z = planes passing through the WD center. The initial temperature 
was taken as T = lO^K everywhere inside the WD. The initial central density was taken as pc = 
2 X 10^ g/cc. Starting from the central pressure P{pc), the equations of hydrostatic equilibrium, 
dP/dr = —GMp/r^ and dM/dr = Anpr'^, were integrated outward until P = was reached. 
The resulting WD configuration was interpolated onto a three-dimensional mesh. Nuclear kinetic 
variables for the unburned matter and / = were defined everywhere. To initiate burning in 
the center, / = 1 was then set inside a small sphere of radius rj = 3 x 10*^ cm and mass Mj ~ 
8 x 10~^ M\YD- For a while, burning was not allowed inside the sphere in order not to impulsively 
disturb the hydrostatic equilibrium. This procedure smoothly triggered the flame propagation 
process just outside the sphere. Due to the finite numerical resolution, the shape of the initial 
flame bubble was not perfectly spherical but contained a number of perturbations. A quadrupole 
mode was the largest mode of perturbations consistent with the symmetries of the problem. Smaller 
perturbations were mostly burned out by the flame, but the quadrupole mode partially survived 
and gave rise to the development of the RT instability. 

3. Results 

3.1. Hydrostatic equilibrium 

Before calculating flame propagation, tests were made to ensure that without burning, the WD 
remains hydrostatic and keeps its symmetry. Since a Chandrasekhar-mass WD is close to a collapse 
threshold, its equilibrium is very sensitive to the discretization errors. Thus a rather fine mesh is 
reqTiircd near the WD center to maintain equilibrium. The initial mesh was constructed with fine 
cells at level / = Imax = 10 for radii r < O.ARwo, at I = Imax — 1 for 0.4:Rwo < r < O.SRwD, at 
I = Imax — 2 for O.SRwD <r< l.lRwD^ and coarser cells at larger r. Computations of hydrostatic 
equilibrium with Ar(Zrreaa;) = 5 x 10^ cm showed numerical noise < 3 x lO^cm/s everywhere except 
in the very outer layers of the WD where the density changed from cell to cell by more than an 
order of magnitude. Near the center, the noise was less than the fiuid velocity generated by the 
laminar flame, ~ Si{pq — pi)/po — 2 x 10^ cm/s. The total kinetic energy of the WD remained 
Ek < 10~^ Eth, where Efh is the total internal energy, during more than two sound crossing times 
of the WD (~ 0.9 s of integration). Deviations from spherical symmetry were negligible. 

3.2. Deflagration of a 0.5C + 0.5O WD 

Deflagration of a CO-WD with initial composition Xc = 0.5 by mass {Yc = 1/24), central 
density pc = 2xl0^g/cc, and angle-independent small scale burning si = 1, S2 = and a = 1 in (7), 
serves as the baseline case. The computational domain for this simulation had size L = 5.35 x 10^ 
cm with levels of refinement Imax = 10, Imin = 7, and with maximum and minimum resolution 
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Ar{lmax) = 5.22 X 10^ and Ar{lmin) = 6.18 x 10^ cm, respectively. The simulation was carried 
out for ~ 1.8s of physical time using Courant number 0.7, required 9,000 timesteps, and took 
~ 6,000 CPU hours. The initial fixed mesh described in Section 3.1 was used for approximately 
Is. After that, the mesh refinement was turned on. The number of cells used in the simulation 
grew from ~ 5, 000, 000 in the beginning to ~ 30, 000, 000 at the end of the simulation. 

Figure 1 presents the evolution with time of the released nuclear energy E^, kinetic energy Ek, 
binding energy Etot, and burned mass Mf, inside the WD. The figure shows that the deflagration 
resulted in a SNIa explosion. At 1.8s, about 1.3 x 10^^ ergs has already been released and Etot — 
8 X 10^'^ergs is positive (the WD is unbound). At this moment, the WD radius is Rwd — 4.7 x 10^ 
cm, pc — 7 X lO^g/cc and the expansion velocity Ugxp — 1-0 x 10^ cm/s. Burning will probably 
release a few more units of 10^° ergs until central density reaches ~ lO^g/cc and burning quenches 
completely. 

Figure 2 shows the time evolution of the flame surface. Initial perturbations lead to the 
formation of rising plumes of burned material. In the beginning, the flame surface and the energy 
release are small, and the flame growth takes place in essentially hydrostatic conditions. At i < Is, 
the kinetic energy remains less that ~ 10^^ ergs or ~ 0.01% of the binding energy of the WD 
(Figure lb). The flame surface increases signiflcantly and becomes more wrinkled at later times. 
The energy release rate increases, and the WD begins to expand. 

Figures 3 and 4 show the flame variable / and the distribution of radial velocity K at 1.44s 
in three orthogonal planes, x = 3.2 x lO^cm, y = 3.2 x lO^cm, and z = 3.2 x lO^cm, offset from 
the WD center. At this time the flame has approximately reached the half-radius of the WD. The 
figures illustrate several important features of burning. Outer parts of the flame consist of big 
plumes rising through the WD. The absolute velocity of these plumes is ~ 6 x lO^cm/s or ~ a^. 
This is about two times larger than the overall expansion velocity of matter, Uexp — 3 X lO^cm/s. 
The plume velocity relative to the expanding material, ~ O.Sog, is subsonic. Figure 4 also shows 
regions of unburned matter sinking to the center. There is a complex flow pattern close to the 
WD center where fresh fuel is burned by the flame and different parts of the flame collide with 
each other and annihilate. This is similar to turbulent burning in a vertical column in a uniform 
gravitational fleld (K95). In that case, the leading part of the flame is also dominated by large 
scale bubbles whose motions determine the amount of fuel passing through the flame. The fuel was 
burned inside the flame brush which is dominated by smaller-scale structures. Since the scale of 
the largest plumes in a column is limited by the column width, the flame is able to reach a steady 
state on all scales. In a spherical star, however, more space becomes available for the plumes as 
the flame moves to larger radii, and the rate of burning must increase with time. 

Figure 5 shows the effective turbulent flame speed St which was estimated as the speed of an 
equivalent spherical flame burning matter with the same rate, 

- ' (15) 



Atrprf dt 

where r/ and p/ are the average radius and the density ahead of the flame brush. In the simulation, 
a small-scale flame speed S determined by (7) was {1 — 2) x lO^cm/s, depending on the radius, 
local density, and time. Figures 1 and 5 show that burning rate increased at t < 1.4 — 1.5s and 
noticeably exceeded the value of S due to the growth of the flame surface. However, St reached 
maximum at ~ 1.5s and started to decline at later times. The WD becomes unbound soon after 
maximum St is reached (see Fig. 1). 

Figures 6 and 7 show the flame and radial velocity at the end of the simulation, 1.8s. The 
leading part of the flame is still dominated by large plumes, whereas inner parts are wrinkled 
and contain more complex structures. The outer edge of the flame is located at ~ Q.^Rwd- 
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The absolute velocity of the plumes has increased and reached lO^cm/s. However, the average 
expansion velocity of the WD has increased significantly and is also Uexp — lO^cm/s. The velocity 
of the flame bubbles is essentially zero relative to the expanding matter, and they have practically 
stopped rising. The expansion completely froze the RT-instability on large scales. Generation of 
new flame surfaces on large scales ceased. Freezing through expansion is the reason for decrease in 
burning rate at t > 1.5s. 

There was an earlier three-dimensional computation of a deflagration explosion in a CO- WD 
(K95) that used a nonuniform expanding 185^ grid, assumed instantaneous transition of burned 
matter to the state of NSE, and the same prescription for S as in this paper. This computation 
was terminated at ~ 1.7s. At ~ 1.7s, the simulation showed plumes of burned material rising at 
c± 5 X lO^cm/s. Approximately 35% of the matter was burned, in agreement with this paper (Note 
that Mf, shown in Figure 16 of K95 must be multiplied by 8, as it erroneously shows the amount 
of burned matter as per one octant of the star) . In the earlier calculation, the WD is still slightly 
bound. 

To illustrate mesh refinement, Figure 8 shows the mesh in the X-Y plane z = 3.2 x lO^cm 
at 1.8s. Figure 9 shows the mesh, density and velocity field in a small region in the X-Z plane 
y = 3.2 X lO^cm with coordinates x = — 1.2 x lO^cm, z = — 9.4 x 10 ''cm. Despite the complexity 
of the fiame surface and velocity field inside the fiame brush, the overall expansion of the WD is 
remarkably one-dimensional. This is due to the subsonic nature of burning which allows pressure 
to equilibrate inside the WD. Small deviations of the WD surface from spherical symmetry are 
only noticeable at late times when the flame gets very close to it. 

3.3. Importance of small-scale burning 

Figure 10 shows the time evolution of E^, Ek, E^ot, and Mf, for the simulation of a 0.5C-I-0.5O 
WD with Ssub = in (6). All other parameters were the same as in the baseline case presented 
in Section 3.2. Flame development was similar in the beginning. However, in this case, burning 
slowed down and then practically stopped because the expansion decreased 5; by ~ 3 orders of 
magnitude. There was no explosion in this case and the WD remained bound. This shows that 
some amount of burning on unresolved scales must be incorporated. Although S ~ lO^cm/s in the 
baseline case was too small to produce an explosion by itself, it is necessary to allow nuclear energy 
release at a flame front whose area increases due to turbulence. A finite S is also necessary to mimic 
the reduction in flame surface caused by merging and burnout of small-scale structures. This is 
crucial in order to reproduce numerically a self-similar regime of turbulent burning in which the 
burning rate becomes independent of the details of flame behavior on small scales (K95). Figure 
11 compares the flame surface at ~ 1.3s for both cases in the X — Y plane passing through the 
center of the WD. This figure illustrates the tendency of the flame to create more surface when S is 
small and to create less surface when S increases. Simulations with more resolution and different 
prescriptions for S are required to determine if we are indeed beginning to see the self-similar 
behavior of the fiame. 

3.4. Nucleosynthesis 

Figure 12 shows the distribution inside the WD of Fe-group, Si-group (intermediate mass 
elements), and unburned CO matter at 1.8s. The location of Si-group elements is determined 
using the criterion 5nse < 0.95 (equation (10)). The figure shows that unburned CO penetrates 
deep inside the fiame brush and is present in significant amounts near the WD center. Incomplete 
burning and the formation of intermediate mass elements occur at densities less than c± lO^g/cc 
and take place at radii greater than c± O.SRwd- Closer to the center, the density is still too high. 
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and burning produces of Fe-group nuclei. The production of S-group nuclei will move inwards as 
the WD expands further. Pockets of completely unburned CO are hkely to survive in central parts. 
Before burning quenches, one should also expect the formation of large amounts of magnesium in 
the central parts of the WD. In one-dimensional models magnesium only forms immediately before 
quenching of incomplete burning in the outer parts of the WD. here, outermost parts of the WD 
are expected to stay completely unburned since the penetration of the flame into the these layers 
has stopped (Section 3.2). 

Figure 13 shows the distribution at 1.8s of the neutron excess of matter r/ = 1 = 2Ye. The max- 
imum neutron excess, rj 0.013, is much smaller than ~ 0.1 typically obtained in one-dimensional 
SNIa simulations (Nomoto et al, 1984; Woosley & Weaver, 1986; Khokhlov, 1991). This is mostly 
due to a reduction of the theoretical electron capture rates (Section 2.4). In addition, burned mat- 
ter spends less time near the center due to RT-mixing. Figure 13 shows that the high-r? products 
are present at radii up to ~ O.SRwd- 

3.5. Sensitivity to initial conditions 

Figure 14 shows the results of a numerical experiment in which the location of the initial flame 
bubble was offset from the WD center by | of its radius (~ lO^cm). The computational domain 
was doubled, and a simplified flame propagation scheme with constant qf = qnse — 9n(0) equal to 
the instantaneous NSE energy release was used. Equations (8)-(ll) were turned off to save time. A 
30% deviation of the ignition cite from the center was enough to cause a highly asymmetric flame. 

In another numerical experiment the initial carbon mass fraction Xc was decreased from 
Xc = 0.5 (Case B) to Xq = 0.2. This resulted in less energy release and less buoyancy of burned 
products. Comparison between the Xq = 0.5 and Xq = 0.2 runs (Figure 15) shows a significant 
delay and slower increase of burning. 

4. Discussion and conclusions 

This paper described a physical model and an adaptive mesh refinement numerical algorithm 
for three-dimensional modeling of Type la supernovae. It also presented the first results of three- 
dimensional simulations using this model of a deflagration in a nonrotating, Chandrasekhar-mass, 
carbon-oxygen white dwarfs. 

The simulation of a 0.5C-I-0.5O CO- WD with central density p = 2x lO^g/cc was carried out 
for ~ 1.8s. The leading edge of the flame was dominated by large-scale, rising plumes of burned 
matter. Between these plumes, unburned CO was sinking towards the WD center. Fine-structured 
flame surface developed closer to the center below the plumes. The simulation lasted long enough 
to demonstrate the effect of freezing of large-scale turbulence due to expansion predicted in (K95) . 
The effective turbulent burning speed St increased exponentially due to the growth of the RT- 
instability, reached maximum value St — 2 x 10^ cm/s at ~ 1.5s, and then started to decline. The 
expansion velocity of the WD reached ~ 10, OOOkm/s by the end of the simulation, and the motion 
of the buoyant plumes relative to unburned matter practically stopped. 

The deflagration resulted in a rather healthy explosion. About 60% of the WD was burned 
by the end of the simulation, and this released ~ 1.3 x lO^^ergs of nuclear energy. The total 
energy, Etot — 8 x lO^^ergs, was positive and the WD was unbound. By 1.8s, the WD radius was 
c± 4.7 X lO^cm, and burning still continued in its central parts. It should be expected that a few 
more units of lO^'^ergs would be released before burning would be quenched completely, so that 
the total kinetic energy at infinity should be in the range (1 — 1.3) x lO^^ergs. This is ~ 50 — 80% 
of the kinetic energy required to produce a typical SNIa explosion. The total burned mass in this 
simulation is expected to be around ~ 1M0. Detailed nucleosynthesis postprocessing is required 
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to predict the amount of ^^^Ni synthesized during the explosion. By analogy with one-dimensional 
simulations, we expect up to ^ 50% of the burned mass, or ~ 0.5Mq, to be ^^Ni. This is nearly 
enough to power a SNIa light curve. 

There are several features of the computed explosion, however, that lead us to conclude that 
it is an incomplete model of a SNIa. One distinct featTirc is the presence of a massive, c± OAMq, 
outer layer of unburned CO surrounding burned matter. The presence of such a layer was a typical 
feature of all one-dimensional deflagration models with constant St or with St decreasing at the end 
of the explosion (Woosley Sz Weaver, Khokhlov 1991). The exception is the W7 deflagration model 
(Nomoto et al., 1976) which created intermediate mass elements at velocities up to c:; 14, OOOkm/s 
by maintaining a high deflagration speed and simultaneously avoiding rapid expansion of the outer 
layers. The simulations presented here show the decrease of the burning rate in the outer layers and 
are not consistent with the W7-like behavior. The presence of the unburned outer layer significantly 
limits the maximum expansion velocities of the intermediate mass elements, and this presents a 
difficulty for the spectra of SNIa at maximum and before maximum light. The observations and 
spectral analyses indicate that the line-forming region for the intermediate mass elements extends 
up to c::! 15, 000 and in some cases up to ~ 30, OOOkm/s in velocity space (Branch, 1981; Benetti et 
al. ,1991; Wells et al., 1994; Hoflich, 1995; Hoflich & Khokhlov, 1996; Lentz et al. 2000). 

The simulations also show that the formation of intermediate mass elements in a three- 
dimensional deflagration explosion occurs throughout almost the entire WD, and that some inter- 
mediate mass elements may form very close to the WD center. This presents another difficulty for 
spectral modeling. The spectra of normal bright SNIa show that the minimum velocity of the line- 
forming region of Si is ~ 8, 000-11, OOOkm/s (Hoflich k Khokhlov, 1996), Ca, ~ 4, 000-6, OOOkm/s 
(Fisher et al., 1997), Mg, ~ 13,000- 14, OOOkm/s (Meikle et al., 1996; Hoflich, 1997; Wheeler et 
al., 1998; Gerardy et al., 2000). There is also an indication that the minimum velocity of carbon 
in SN1990N (CII at 6580 A) is larger than 26,000 km/s (Fisher et al. 1997), which means that 
practically all of the material in this supernova has been burned. According to the simulations, 
unburned CO is also likely to remain in the central regions of the deflagration supernova where 
it will coexist with ^^Ni. Whether this can lead to the excitation of oxygen lines in the late- 
time spectra when the positron energy deposition becomes nonlocal (Milne et al. 1999)is an open 
question which requires further investigation (Hoflich, private communication). 

The formation of a massive outer layer of unburned material is a large-scale effect caused by 
the competition between the of buoyancy of large-scales structures (rising plumes) and the global 
expansion of the WD. In this competition, buoyancy eventually loses. During the explosion, the 
expansion velocity of the star increases gradually, and at 1.8, outer layers have velocities comparable 
or even exceeding the velocity of the large-scale plumes of burned products. Even if burning on 
small scales is underestimated in the simulation, increasing the energy released by small scales 
would increase the expansion velocity, but would hardly affect the buoyancy of the plumes. As the 
WD continues to expand into the vacuum, the rarefaction will accelerate the outer layers and this 
will further increase their velocity relative to the burned material. 

The difficulties associated with the deflagration explosion can be avoided if burning turns 
into a detonation during the explosion (delayed detonation). Detonation would easily overcome 
the expansion and would incinerate the outer unburned layers of CO. This was one of the initial 
reasons for the introduction of the delayed detonation model of SNIa. It is also clear that the 
delayed detonation will help with the potential problems associated with three-dimensional effects 
by wiping out composition inhomogeneities in the inner layers. 

These are the flrst simulations with the model presented in the paper. As such, they show 
new effects but also rais many questions. Computations with more resolution and with different 
subgrid burning models are required to determine if the resolved scales are small enough that 



11 



burning on these scales is in a self-similar regime. In particular, at maximum burning (2:^ 1.5s), 
the simulations show large-scale motions with relative velocities of ~ 5 x lO^cm/s. Thus there 
should be shear instabilities and a strong turbulent cascade on the sides of rising plumes. This 
should be incorporated into the model. The question of the deflagration-to-detonation transition 
(DDT) must also be addressed. A rough estimate can be made as follows. Assume that at 1.5s, 
the turbulent velocity is ~ 10*cm/s on the driving scale L ~ lO^cm. Then, from Kolmogorov 
cascade, the turbulent velocity U\ ~ 10^ A cm/s on a scale A, where A is in centimeters. At 
p c± lO^g/cc, the thickness of the laminar flame is 5i ~ lO^cm, and on the scale A ~ (5/ the 
velocity U\ lO^cm is much greater than Si. Thus, the turbulent velocity on this scale may be 
strong enough to move the flame from the flamelet into the distributed regime, create hot spots 
of nonuniformly preheated unburned material, and thus trigger a detonation (Khokhlov et al., 
1997). To trigger a detonation, the size of a nonuniform region must exceed a critical value which 
depends on density. At p ~ lO^g/cc, the critical size lO'^cm is much less than RwD and the size 
of large scale structures obtained in the simulations. Thus the transition to detonation seems to 
be possible. Significant work requires to test these ideas. This paper shows that high-resolution 
three-dimensional modeling of SNIa is now a reality. 

The calculations presented here indicate that the result of the explosion may be sensitive to 
varying the position of the ignition site and the initial carbon fraction Xq in a WD. Decreasing of 
Xc leads to less energy release and less buoyant products of burning. This results in much slower 
development of the RT-instability, so that the freezing due to expansion may be more effective. 
On the other hand, smaller energy release leads to slower expansion as well. It is most likely that 
decreasing Xc would result in a weaker explosion. Sensitivity to the position of the ignition site 
indicates that the results may also be sensitive to stellar rotation, which may lead to a preferential 
growth or suppression of the initial perturbations along the rotational axis. This makes both Xc 
and rotation potentially good candidates for being the 'hidden parameters' to describe the observed 
diversity of SNIa (Phillips, 1993). 
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Figure captions 

Figure la. - Figure la. - Total (Etot), kinetic {Ek), released nuclear (En) energy and burned mass Mb of 
a |C + WD during a deflagration explosion (Section 3.2). 

Figure lb. - Solid line - kinetic energy (E^) of a |C + WD during a deflagration explosion at earlier 
times (Section 3.2). Dashed line - kinetic energy with burning turned off (Section 3.1). 

Figure 2a. - Time evolution of the flame surface during a deflagration explosion (Section 3.2). 

Figure 2b. - Same as Fig. 2a. Details of the flame structure at 1.8s. 

Figure 3. - Flame variable / (Equation (2)) for a deflagration explosion Section 3.2 at 1.44s shown in 
three orthogonal planes offset by 3.2 x 10''cm from the WD center. 

Figure 4. - Same as Figure 3 but shows the radial velocity. 

Figure 5. - Estimated turbulent flame speed St (Equation (15)) during a deflagration explosion Section 
3.2. 

Figure 6. - Flame variable / (Equation (2)) for a deflagration explosion Section 3.2 at the end of the 
simulation 1.8s shown in three orthogonal planes offset by 3.2 x lO^cm from the WD center. 

Figure 7. - Same as Figure 6 but shows the radial velocity. 

Figure 8. - Mesh and flame variable / in the X-Y plane passing through the WD center for a deflagration 
explosion Section 3.2 at the end of the simulation t = 1.8s. 

Figure 9a. - Density and velocity in the X-Z plane passing through the WD center for a deflagration 
explosion Section 3.2 at the end of the simulation 1.8s. The region x = 1.2 x lO^cm and 
z = — 9.4 X lO^cm is shown. 

Figure 9b. - Mesh in the X-Z plane passing through the WD center for a deflagration explosion Section 
3.2 at the end of the simulation 1.8s. The region x = 1.2 x lO^cm and z = — 9.4 x lO^cm is 
shown. 

Figure 10. - Time evolution of the total {Etot), kinetic {Ek), released nuclear {E^) energy and burned mass 
Mb of a + |0 WD during a deflagration explosion with the laminar flame speed S = Si 
(Section 3.3). 

Figure 11. - Flame surface in the X-Y plane at c± 1.3s for the simulation with (right panel, Section 3.2) 
and without (left panel, Section 3.3) turbulent small-scale burning. 

Figure 12. - Distribution of the Fe-group, Si-group, and CO matter for a deflagration explosion Section 3.2 
at the end of the simulation t = 1.8s shown in three orthogonal planes offset by 3.2 x lO^cm 
from the WD center. 

Figure 13. - Same as Figure 12 but shows the neutron excess r] = 1 — 2Ye. 

Figure 14. - Deflagration starting with the offset ignition point (Section 3.5). Shows the flame variable and 
velocity at c± Is. Black square shows the ignition point. 

Figure 15. - Evolution of the kinetic energy during deflagration explosion of a WD with the initial 0.5C + 
0.5O and 0.2C -I- 0.8O compositions (Section 3.5). 
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